/** 
 * \file projections.c
 *       Filename:  projections.c
 * =====================================================================================
 * 
 *    Description:  Implementation of fonctions related to cartographic projection
 * 
 *        Version:  1.0
 *        Created:  11/11/2009 16:08:09
 *       Revision:  none
 *       Compiler:  gcc
 * 
 *         Author:  François Hissel (fh), francois.hissel@m4x.org
 *        Company:  
 * 
 * =====================================================================================
 */

#include	<math.h>
#include	"projections.h"

double convert_angle(double deg,double min,double sec) {
	return (deg+min/60+sec/3600)*M_PI/180;
}

void _molodensky(double u,double v,double *x,double *y,double *params) {
	double s1=sin(v);
	double s2=s1*s1;
	double e2=2*params[3]-params[3]*params[3];
	double t=sqrt(1-e2*s2);
	double n=params[2]/t;
	double rho=params[2]*(1-e2)/(t*t*t);
	*x=u+(-params[4]*sin(u)+params[5]*cos(u))/(n*cos(u));
	*y=v+(-params[4]*sin(v)*cos(u)-params[5]*sin(v)*sin(u)+params[6]*cos(v)+sin(2*v)*(params[2]*(params[3]-params[1])+params[3]*(params[2]-params[0])))/rho;
}

void molodensky(double u,double v,double *x,double *y,char code) {
	double params[7];
	switch (code) {
		case 0:
			params[0]=6378249.2;
			params[1]=1-6356515.0/6378249.2;
			params[2]=6378137.0;
			params[3]=1/298.257223563;
			params[4]=-168;
			params[5]=-60;
			params[6]=320;
			break;
		case 1:
			params[0]=6378137.0;
			params[1]=1/298.257223563;
			params[2]=6378249.2;
			params[3]=1-6356515.0/6378249.2;
			params[4]=168;
			params[5]=60;
			params[6]=-320;
			break;
		default:
			return;
	}
	_molodensky(u,v,x,y,params);
}

void _lambert(double a,double b,double *x,double *y,double *params) {
	double e=sqrt((params[0]*params[0]-params[1]*params[1])/(params[0]*params[0]));
	double l=log((1+sin(b))/(1-sin(b)))/2-e/4*log((1+sin(b)*e)/(1-sin(b)*e));
	double r=params[4]*exp(-params[3]*l);
	double g=params[3]*(a-params[2]);
//	double n=log(cos(params[2])/cos(params[3]))/log(tan(M_PI/4+params[3]/2)/tan(M_PI/4+params[2]/2));
//	double f=(cos(params[2])*pow(tan(M_PI/4+params[2]/2),n))/n;
//	double rho=f/pow(tan(M_PI/4+b/2),n);
//	double rho0=f/pow(tan(M_PI/4+params[1]/2),n);
	*x=params[5]+r*sin(g);
	*y=params[6]-r*cos(g);
}

void lambert(double a,double b,double *x,double *y,char code) {
	double params[7];
	params[0]=6378249.2;
	params[1]=6356515.0;
	params[2]=convert_angle(2,20,14.025);
	switch (code) {
		case 0:
			params[0]=6378137.0;
			params[1]=params[0]*(1-1/298.257222101);
			params[2]=convert_angle(3,0,0);
			params[3]=0.7256077650;
			params[4]=11754255.426;
			params[5]=700000;
			params[6]=12655612.050;
			break;
		case 1:
			params[3]=0.7604059656;
			params[4]=11603796.98;
			params[5]=600000;
			params[6]=5657616.674;
			break;
		case 2:
			params[3]=0.7289686274;
			params[4]=11745793.39;
			params[5]=600000;
			params[6]=6199695.768;
			break;
		case 3:
			params[3]=0.6959127966;
			params[4]=11947992.52;
			params[5]=600000;
			params[6]=6791905.085;
			break;
		case 4:
			params[3]=0.6712679322;
			params[4]=12136281.99;
			params[5]=234.358;
			params[6]=7239161.542;
			break;
		case 5:
			params[3]=0.7289686274;
			params[4]=11745793.39;
			params[5]=600000;
			params[6]=8199695.768;
			break;
		default:
			return;
	}
	_lambert(a,b,x,y,params);
}

void projection(double a,double b,double *x,double *y) {
	double u,v;
	molodensky(a,b,&u,&v,1);
	lambert(u,v,x,y,5);
}
